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The cavity method is a well established technique for solving classical spin models on sparse 
random graphs (mean-field models with finite connectivity). Laumann, Scardicchio and Sondhi 
arXiv : 0706 . 4391 proposed recently an extension of this method to quantum spin-1/2 models in a 
transverse field, using a discretized Suzuki- Trotter imaginary time formalism. Here we show how to 
take analytically the continuous imaginary time limit. Our main technical contribution is an explicit 
procedure to generate the spin trajectories in a path integral representation of the imaginary time 
dynamics. As a side result we also show how this procedure can be used in simple heat-bath like 
Monte Carlo simulations of generic quantum spin models. The replica symmetric continuous time 
quantum cavity method is formulated for a wide class of models, and applied as a simple example 
-^ , on the Bethe lattice ferromagnet in a transverse field. The results of the methods are confronted 

with various approximation schemes in this particular case. On this system we performed quantum 
Monte Carlo simulations that confirm the exactness of the cavity method in the thermodynamic 
limit. 



I. INTRODUCTION 
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Mean-field approximations are often useful first steps to unveil the physical content of realistic models. This is all 
the more true when exact solutions are probably impossible to obtain in the finite-dimensional setting, in particular 
p^ \ when quenched disorder and/or quantum effects have to be taken into account, as for instance in the case of Anderson 
^ ■ localization!. Another example is dynamical mean- field theory^, that has been a very fertile approach to the problem 
of strongly correlated fermions. ft can be sometimes preferable to study mean-field theories not by making an 
approximation to a finite-dimensional model, but rather by formulating a model which is mean-field by nature, this 
allowing in particular to state the results in a mathematically clearer way. The simplest of such examples is the 
Curie- Weiss model of ferromagnetism, in which N classical Ising spins all interact attractively which each other, with 
a coupling constant scaling inversely with the size of the system to ensure a well-defined thermodynamic limit. The 
equivalent for quenched disordered systems is the Sherrington- Kirkpatrick model^ of a spin glass, where again all 
spins interact weakly with other, yet with coupling constants of random signs. 

The mean-field character of the above mentioned models arises from their infinite connectivity (in the thermody- 
namic limit). There exists however another class of models, which are still mean- field yet keep a finite connectivity, 
each of the degrees of freedom they possess interacting only with a finite (with respect to N) number of neighbors. For 
ferromagnetic models they can be obtained by the Cayley tree construction, where one draws an infinite regular tree 
and studies the magnetization of the root sites. Cayley tree models have however pathological surface effects, and the 
theory of finitely-connected mean-field frustrated systems is better defined on random graphs 5,6 , of fixed or fluctuating 
connectivity. Classical models of spins on such random structures have been the subject of extensive study in the last 
decade. These works were motivated on the one hand by their somehow more physically realistic features, namely the 
finite connectivity, and also because of their strong relationship with issues originating from computer science, namely 
the understanding of phase transitions in random constraint satisfaction problems^^. Finite-connectivity models 
are technically much more involved than their fully-connected counterparts. The replica method^ that has been first 
developed to solve the Sherrington-Kirkpatrick model becomes less practical in this setting^, and the alternative 
cavity method turned out to be more useful 6 . 

The interplay between quenched disorder and quantum fluctuations can lead to a very rich phenomenology, and 
in particular the properties of the glass phase found at low temperatures in classical models can be qualitatively 
modified when a transverse field acts on the system^. More generally the issue of the nature of the quantum phase 



transitions at zero temperature^ in presence of disorder is a very rich one. In the context of mean-field theory this 
point has been mainly studied in fully-connected model a 14 ' 15 ' 16 ' 17 ' 18 ' 19 ' 20 , with a few exceptions that appeared in 
the last yea r 21 i 22 ' 23 i 24 . In a very interesting contribution^! Laumann, Scardicchio and Sondhi made a first step in 
extending the cavity method to quantum spin models in a transverse field, and in this paper we shall develop further 
this idea by solving a discretization problem which plagued their proposal. Let us also mention here the work of Knysh 
and Smelyanskiy 23 , who developed a similar approach in the framework of the so-called static approximatio n 14 ! 16 ' 18 . 
The motivations for this line of work is two-fold. On the one hand one can expect an even richer physical behavior 
of finitely-connected quantum models with respect to the fully-connected ones. The possible fluctuations in the local 
geometry, and some notion of distance which was absent in fully-connected models opens the way to a more complex 
phenomenology. On the other hand one can aim at a better understanding of some issues of quantum computing, and 
in particular on the use of quantum annealing (or adiabatic algorithm) 25 ' 26 ' 27 to solve random constraint satisfaction 
problems. These quantum algorithms do indeed rely on the application of a transverse field on spin models that have 
been extensively studied at the classical level with the cavity method. Computing the location and nature of the 
phase transitions^ 8 - encountered along the annealing path (as the transverse field is progressively turned off) might 
give some informations on the behaviour of these quantum algorithms themselves. 

The remaining of the article is organized as follows. In Sec. |TT] we recall the Suzuki- Trotter approach to spin-1/2 
models in a transverse field, and develop our main technical contribution in Sec. Ill Bl where we show how to actually 
build the spin trajectories of the path integral representation of the imaginary time evolution operator. Section [Ml is 
then devoted to the study of a very simple example of finitely-connected quantum model, namely the Bethe lattice 
ferromagnet. We first explain the continuous time quantum cavity treatment of this model, before presenting the 
results of the method and confronting them with some approximate approaches. In Sec. IIVI we present numerical 
results of Monte Carlo simulations we performed for this model, and show how the computations of Sec. Ill Bl can be 
turned in a simple and versatile quantum Monte Carlo method. The generic formalism of the quantum cavity method 
is developed in Section [V] we hope this order of presentation, and the inclusion of a fully worked-out example before 
the general case, will ease the reading of this work. We finally draw our conclusions and put forward perspectives for 
future work in Sec. IVII Some technical details are deferred to a series of Appendices. 



II. PATH INTEGRAL REPRESENTATION FOR SPIN-1/2 MODELS 

A. Spin models in a transverse field: Suzuki-Trotter formalism 

Let us consider the Hilbert space spanned by the orthonormal basis of 2 N kets |er), where a = (ci, . . . , &n) denotes 
a configuration of N Ising spins, er^ = ±1. This space can be viewed as the tensorial product of N spins 1/2, with 
operators a\ and erf , whose action on the base vectors is defined by 

of Iff) = *i|2) . (1) 

cf|cr) = |<7l,. ..,CTj_i,— <Ti,(Ti+l,...,<Ttf) ■ 

From a classical energy function of N Ising spins, E{a\, . . . , <tn), one can construct an operator E — E(af, . . . , a z N ), 
diagonal in the {|ct}} basis. The Hamiltonian operators investigated in this paper are obtained from such a classical 
energy by the addition of a transverse field, 

N 

H = E-B^o* , withB>0. (2) 

Our goal is then to compute the quantum statistical mechanics properties at inverse temperature (3, i.e. the partition 
function Z and the average of observables (operators) O, defined by 

Tr (6 e-P") 
Z = Tr (e- pH ) , (6) = \ -^- . (3) 
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A well-known way of tackling such problems is to transform them into an extended Ising model by using the Suzuki- 
Trotter formula"*, as summarized in the following lines : 



Z = Tr((e-%* + % B ^°<) NB ) (4) 



= lim Trf ( e --^ S e^ B ^= iaf ) 

JV s ^oo y V / 

= lim Y TT (a a \e-^ S e^ B ^°'\a a+1 ) 

iV s — *oo *- — ' -*--*- 
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In the two last lines g_ Ns+1 = a . For a finite value of the number of Suzuki- Trotter "slices" N s , the problem 
has thus become one of N X N s Ising spins, each of the at being promoted to a ring (aj, . . . , a i B ) with nearest 
neighbor ferromagnetic interactions along the "discrete imaginary time" a axis (with periodic boundary conditions). 
The original interactions E acts indentically and independently on each of the configurations q_ a . For notational 
convenience we shall use bold symbols for quantities that depend on the slice a, for instance <Ji — {a}, . . . ,a i B ) is the 
configuration of the ring of Ising spins at site i and g_ = (er 1 , . . . ,a N ") is the full configuration of the N x N s spins. 
We can thus introduce a probability measure on the N x N s Ising spins, 



/*(2) = ^e-^f]^), (5) 

AT 
er_ i=l 

such that the normalization constant Z^ B reduces to the partition function Z in the N s — * oo limit (in the following we 
shall sometimes keep implicit the dependence on N s ). To write in a compact way this last equation we have defined 

1 "• 
%)^E%°)- (6) 

s a— 1 

the average of independent copies of the classical energy on the various slices, and 

<°) = n<a Q |e^V +1 ) = fl ( cosh (ip) S^+i+swh (?p) ^,_ CT « +1 ) , (7) 

a=l a=l \ \ s / \ s / / 

the ferromagnetic interaction along the imaginary time axis induced by the transverse field (we use a Ns+1 —a 1 ). One 
can easily show that the average value of observables can be obtained in this formalism as 

(O) = ]T m (a: — rw — — " > ( 8 ) 

where the slice number a is here arbitrary, thanks to the cyclic invariance around the discrete imaginary time axis. 
This can be simplified further for observables O diagonal in the {\a}} basis, i.e. that can be written as 0(af, . . . , cr%j)' 



, N. 
(O) = £ v(<L)0(a a ) = J2 »{<L) W E °^) ■ ( 9 ) 



a=l 



A non-diagonal observable we shall study in the following is the transverse magnetization (erf) (written here for an 
arbitrary site i), which can be computed as 

(of) = V^)f E ^ '" T - } =E^)F E ftanh (P* 



B. The continuous imaginary time limit 

To recover the truly quantum properties of the model one has to perform the limit N s — » oo. The basic degrees 
of freedom cr^ which were the configurations of a ring of Ising spins [a] , . . . , a i ") then becomes piecewise constant 
functions cr^t) £ { — 1,1} of an imaginary time parameter t, the discrete coordinate a S [1, JV S ] being mapped to 
t E [0,(3] with the correspondence t = j3a/N s . In this limit the sum over er in the expression ([5]) of the partition 
function is naturally interpreted as a path-integral. The discreteness of the spin degrees of freedom actually make 
such a path-integral representation— easier to formulate than Feynman path-integrals for continuous coordinates^, 
and can be given a rigorous mathematical conten t 32 i 33 ' 34 i 35 . Note that these continuous time trajectories can be 
easily represented in the memory of a computer, as the trajectory of site i is fully specified by <7j(£ = 0) and the 
times at which the spin flips. Actually numerous continuous time quantum Monte Carlo algorithms do exist, see for 
instanco 36 ' 37 i 38 ' 39 . 

The rest of the paper will crucially rely on the procedure developed in Sec. Ill Bl2l and HT B 31 Though it will also be 
useful for analytical purposes, it is more intuitively motivated by the following simulational consideration. Maybe the 
simplest way to ensure the detailed balance condition in a Monte Carlo simulation which aims at sampling an arbitrary 
measure /u(<r) is to perform transitions from the current configuration a to a configuration obtained by replacing the 
value of a randomly chosen degree of freedom Ui, by a random value drawn from the measure conditioned on all other 
degrees of freedom. This procedure is known in classical simulations as the heat-bath, or Glauber algorithm. Its 
equivalent in quantum simulations consists in drawing a new configuration of the ring <Tj, or of the trajectory <Ji{t) 
in the continuous imaginary time, according to the equilibrium measure induced by the spin trajectories of all other 
sites. A moment of thought reveals that this boils down to study the evolution of a single spin-1/2 in the presence of 
a constant transverse field and a piecewise constant longitudinal field, the latter being the effective field induced by 
the rest of the system on <Ji. This is precisely the issue we shall tackle in Sec. IIIB2l and lHB3| after having recalled 
in Sec. Ill B li the well-established path-integral representation of a spin-1/2. 

1. The path integral representation of a single spin in constant fields 

Let us define the propagator for the evolution during an interval of imaginary time A of a spin in constant transverse 
and longitudinal fields (B and h respectively): 

W(a — a',h,\) = (o-\e x(ha * +Ba ^\cT') . (10) 

The diagonalization of the order 2 matrix ha z + Ba x easily leads to 



[ cosh(AV^ 2 + h 2 ) + a ,J tl .. smhjWB 2 + h 2 ) if a = a 1 



if a = -a' 



The path-integral representation of this propagator read: 



30,38 
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W{<t -Kr,h,\) = Y^ B2n / di i / di 2 ■ • • / dt 2n exp[ah(2ti - It 



2 + 2t 2n + A)] , (12) 



°° P P P 

W(a -» -a, h,X) = V B 2n+1 / dii / di 2 . . . / di 2rl+ i exp[ ( j/i(2i 1 - 2t 2 + • • • + 2t 2 „ +1 - A)] . (13) 

Each term of these expressions corresponds to a spin trajectory that changes value at times t\ < t% < • • • ; it is 
weighted by a factor B raised to the number of such discontinuities, and by exp[/i J Q a(t)]; a spin trajectory with 
identical (resp. opposite) initial and final value has to jump an even (resp. odd) number of times. There are two 
ways to convince oneself of the correctness of this result. Applying the Suzuki- Trotter formalism to this single spin 
problem leads to such a weight in the N s —*oo limit^. Alternatively one can notice that Eqs. (TT2^) . (|13p coincide with 
(fTTjl at A = and that they obey the same set of first order linear differential equations, 

d 

— W(a -> a' ', h, A) = cr'h W(<J -> a', h,\)+B W{a -> -a> , h, A) , (14) 

which implies that they coincide for all values of A. 
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FIG. 1: A pictorial representation of Eqs. (|16[) . (|17l) . 

2. Generating trajectories for a constant longitudinal field 

The above expressions (|12[) . (|13p can be interpreted as the normalizing constants of probability measures on the set 
of piecewise constant functions from t e [0, A] to { — 1, +1}, conditioned on their initial (a(t = 0)) and final (a(t = A)) 
values. More explicitly, for instance for a(t = 0) = a(t = A) = a, the probability of a trajectory with In flips at times 
in the infinitesimal intervals [tj, tj + dtj], with t\ < ■ • • < t^m is defined to be 

B 2n e ah(2t 1 -2t 2 + --2t 2n+ X) ^ _ _ ^ _ (15) 



W(a^KT,h,X) 

Our goal is now to construct a procedure for actually sampling from these probability measures, that is constructing 
spin trajectories according to these weights. We shall do this by exploiting the following two identities: 

W(a -f a, h, A) = e ahx + B f du e ahu W(-a -» a, h, A - u) , (16) 

Jo 

W(a-*-cr,h,X) = B [ du e ahu W(-o -v -a, h, A - u) . (17) 

Jo 

The path-integral interpretation of these relations, more easily conveyed by the drawing of Fig. [TJ is as follows. For 
the first one, it means that a spin trajectory starting and ending at the same value cr(0) = c(A) = a is either constant 
on the whole time interval, or made of a constant part upto time u, followed by a jump to — a and a second part 
of the trajectory representative of W(—a —* a,h,X — u). Similarly the second one expresses the necessity for a 
trajectory from a to —a to have at least one discontinuity at a given time u, followed by a trajectory accounting for 
W(—a — > — a, h,X — u). These equalities can be proven either from the path-integral representation of (p~2|) . (p~3|) . or 
from the explicit expressions of W given in Eq. (jTTJ). As a consequence of (fT6|) . (fl7|) one obtains the following recursive 
procedure to draw a spin trajectory for a constant longitudinal field on a time interval of length A, constrained to 
cr(0) = a, o-(A) = a' : 

• if a = —a' 

— draw a random variable u € [0, A] with density proportional to e ahu W(—a — ► — a, h,X — u) (see below in 
cq. ([T8|) for some details on how to perform this step) 

— set a(t) — a upto time u 

— call the same procedure to generate a trajectory from — a to — a on the remaining interval of length A — u 

• if a = a' 

— with probability e ahx jWia — ► a, h, A), set a(t) = a on the whole time interval, and exit the procedure 

— otherwise, 

* draw a random variable u £ [0, A] with density proportional to e ahu W(—a — * a, h, A — u) 

* set a(t) = a upto time u 

* call the same procedure to generate a trajectory from —a to a on the remaining interval of length A — u 

In order to draw u £ [0, A] with a density proportional to e ahu W{—a —> —<r, h,X — u), we compute its cumulative 
distribution, 



£ dt e° ht W(~a -* -a, h,X-t) = i _ hu sinh((A - n)VB 2 + fe 2 ) 
/ A dt e" ht W(-a ->-<r t h,\-t) sinh(AV5 2 + h? 



G{u) = jo_^__^ -- ^2 <2 = 1 ^^; < ;^ i ;^ . ( i 8 ) 
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FIG. 2: Definition of the effective field trajectory. 

A simple way to draw u amounts to draw G uniformly at random on [0, 1], and to invert the above expression to 
obtain u(G). One can proceed similarly for the generation with the density proportional to e u W(—a — ► a, h, A — it), 
which involves another cumulative distribution G. 

3. Generating trajectories for a piecewise constant longitudinal field 

In the previous subsection we considered the particular case of a constant longitudinal field h. Let us now address 
the general case of a piecewise constant h — h(t) on the interval of imaginary time [0, j3] , with the following definitions 
illustrated in Fig. O we shall call p the number of times it changes value between t = and t = /?, = t^ < t^ < 
■ ■ ■ < t^ < £(p +1 ) = (3 the times of these changes, Aw = £( 4+1 ) — fM the length of these intervals for i g [0,p], and 
finally h^ the values the field takes in each of these intervals. We shall have to compute the partition function of a 
spin acted on by such a field, 

Z(h)=Tr(f[e xm ^ m ^ +Ba 'y\ , (19) 

and to generate spin trajectories according to the corresponding weights. The computation of the partition function 
can be performed by inserting p + 1 representations of the identity in (fT9"|) , corresponding to the spin values at the 
imaginary times t^ l > where h(t) is discontinuous: 

z ( h ) = Yl Z(a ,...,o- p \h) , 

<7o,...,<T p 
V 

Z(a ,...,a p \h) = Y[W(ai^a i+1 ,h^,X^), (20) 



i=0 



with cr p+ i = (To- Given the trajectory h this computation is easily performed, necessiting only the multiplication 
of the p matrices of order 2 defined in (|10[) . The sampling of the spin trajectory a(t) on the interval [0,/3] is done 
as follows. The values (ao, ...,a p ) of the spin at times (t^ ', ■ ■ ■ ,t^) are generated^ according to the probability 
Z(<jq, . . . ,o~p\h)/Z(h). Then for each interval i £ [0,p] a spin trajectory from Ui to ct^+i is generated according to 
the procedure of Sec. Ill B 2[ the longitudinal field being constantly equal to hS % ' on this interval of time. Finally the 
p + 1 trajectories are concatenated to obtain the full trajectory from t = to t = /3. 

Let us emphasize that the path integral representation of the imaginary time evolution is well known in the 
literatur o 30 i 31 ' 32 i 33 ' 34 i 35 ; however we could not find in previous works such an explicit sampling procedure for generating 
the spin trajectories. Actually, as far as we know, all continuous time quantum Monte Carlo algorithm a 36 i 37 ' 38 i 39 do 
not proceed in an heat-bath way, generating "from scratch" a new spin trajectory conditioned on the local effective 
field, but rather constructing the spin update using the current configuration itself. 
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FIG. 3: A pictorial representation of Eq. (|22l) . 

III. THE QUANTUM BETHE LATTICE FERROMAGNET 

The remainder of this article will be devoted to the study of the simplest of the finitely connected models that can 
be handled by the quantum cavity method, namely the transverse field quantum spin-1/2 ferromagnet on the Bcthc 
lattice (more precisely on a random regular graph). The physical properties of such a model are very intuitive: at 
low temperature and transverse field the model is ferromagnetically ordered, with a positive spontaneous longitudinal 
magnetization. Thermal (increasing T) or quantum (increasing B) fluctuations destroy this order outside a region 
delimitated by a critical line in the (B,T) plane, that ends up in a quantum critical point at zero temperature. 

An even simpler model displaying these features is the (fully connected) quantum Curie- Weiss model, whose solution 
we recall in Appendix lAl Both of them are of a mean-field nature, and should share most of their qualitative properties, 
yet the Bethe lattice model is quantitatively different, and technically more involved because of its finite connectivity. 

A. The quantum cavity method treatment 

The quantum Bethe lattice ferromagnet is defined by the Hamiltonian 

N 

H=-J2vi^-Bj2< > ( 21 ) 

i-j i=l 

where the first sum runs over the edges of a random I + 1-regular graph of N vertices* -. This means that the graph of 
interactions is uniformly chosen among all graphs on N vertices for which each vertex has the same number I + 1 of 
neighbors. These graphs have the good properties to realize finite size Bethe lattices: the set of vertices at distance 70 
smaller than a given cutoff d from an arbitrarily chosen vertex is, with a probability which goes to one when N diverges 
with d held fixed, a regular tree of connectivity I + 1. On the other hand such a graph has no surface, contrarily to 
the usual Cayley tree, and the problem of the boundary conditions for frustrated models (that will be encompassed 
by the generic treatment of Sec. |V| is absent with such a definition. 

The hypotheses of the cavity method are more simply explained assuming first that the graph of interactions is 
actually a finite tree, and that the quantum aspects of the problem have been handled by a finite number N s of 
Suzuki- Trotter slices. In such a case it is easy to solve the model exactly by taking benefit of the natural recursive 
structure of a tree: one breaks the graph of interaction into subtrees that are then glued together. Let us explain 
this with more precise formulae, defining /ij_>j (<Tj) as the probability law of the configuration of the ring <Xj when the 
interaction with its neighbor j has been removed from the graph. If we denote by di the set of vertices neighbors of 
i, the recurrence equations for these distributions are: 

Mi-i(o-i) = — w(<n) n ^M^(Tfc)e /3CT *" Tfc , (22) 

i ~* i k£di\j T k 

where %— >j is a normalization constant, and we introduced for two arbitrary imaginary time dependent quantities 
a ■ b = X] a =i a a b a /N s . A graphical representation of this equation is given in Fig. [3] in the case I — 2. The site i 



8 

has then two neighbors, k\ and k%\ the distributions fj,^ 2 ^i encode the effect on ki t 2 of the part of the tree that does 
not involve i; this is represented by the dashed line in figure [3J In absence of site i the sites k\ and fc 2 are decoupled 
(as we assume the graph to be a tree there are no paths between fci and k^ that do not go across i) and hence 
independent. Therefore the distribution of <Tj in absence of site j is given by (|22[) . A proof of such recursion equations 
and a discussion of the connections with the Bethe-Peierls approximation and the Belief Propagation algorithm can 
be found i n 41 ' 42 . On a given tree this set of equations (one for each directed edge of the graph) has a unique solution, 
that is a set of "messages" fii^>j, that can be efficiently determined by sweeping the edges from the leaves towards 
the center of the tree. From them all the relevant thermodynamic quantities can be computed. For instance the 
probability law of the configuration of the ring <Tj in the complete tree is given in terms of the messages sent by its 
neighbors, 

/i(o-i) = -w(*i) Y[ 5^Wfe-i(ffk)e / *"" ff * , (23) 

with Zi a normalization constant, while the joint law for two neighboring sites i and j reads 

n{<T i ,a j ) = m^ j (a-i)n J -,i((T j )e^' Ti " T: > . (24) 

Zi-j 

Moreover the free-energy per site / can be expressed in terms of the normalization constants of these laws, 

1 1 N 1 

_ / 3 / = _lnZ=-5>* i --5>* i _, ) (25) 

i=\ i-j 

where the second sum runs over the (undirected) edges of the tree. 

The above derivation was exact because we assumed the graph of interactions to be a tree. The scope of the cavity 
method is to extend these results to random graphs which are only locally tree-like, in the precise sense explained 
above. In its simplest version, called replica symmetric for historical reasons, one assumes the existence of a single 
pure state in the configuration space of the random graph model. This implies a decay of correlations at large distance 
in the graph, hence the effect of the long loops neglected in the tree derivation amounts to create a self-consistent 
boundary condition which traduces the absence of a surface in the random graph. As in the present model all sites 
have the same neighborhood (there is no fluctuation neither in the connectivity nor in the intensity of the interaction 
couplings), one has to look for a solution of (f2"2")l where the messages /ij->j are all equal to a single law, that we shall 
denote in the following ?y(o"), which is seen from (j2"2")l to satisfy 



(a) = ^-w(a)r£ V (a')e^A 



r)(a) = -w(a) r^^y™ j ■ ( 26 ) 

The assumption on the unicity of the pure state can fail for two kind of reasons. In ferromagnetic models, as the one 
considered now, there is an ordered phase at low temperature and transverse field in which two pure states coexist 
because of the up/down symmetry of the model. This is not a serious limitation of the method, in the following it 
will be kept understood that an infinitesimal longitudinal field is applied to the system in order to select one of the 
two pure states. In fact the exactness of the cavity method for classical ferromagnetic models on random graphs has 
been recently proven rigorously 4 ^. A much more serious problem, that shall not be discussed in this paper, arises in 
frustrated models, when an exponential number of pure states proliferate at low temperatures; the replica-symmetry 
breaking version of the cavity method is then required to solve the problem. 

Let us first write the solution of the equation (|26| in the classical situation, for B = 0. In such a case the weight 
w forces all spins along the ring to take the same value, hence 

m = '-^^ (ft w) + '-^^ { fi w.) , P7) 



\a— 1 / \ck— 1 



where h is solution of 



h= -arctanh(tanh(/3)tanh(/3/i)) . (28) 

P 

The classical model thus exhibits a continuous paramagnetic to ferromagnetic transition when the inverse temperature 
P crosses its critical value, /3 C = arctanh(l/Z). 
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FIG. 4: Longitudinal (m z , left panel) and transverse (m x , right panel) magnetizations as functions of the transverse field B for 
T — 1.7, 1.0, 0.1. For the lowest temperature we could not obtain reliable data for m z close to the critical field due to critical 
fluctuations. The dashed line is a fit to m z ~ (B c — B) 1 ' 2 done in the interval B € [2.1, 2.2] and serves as a guide to the eye. 

We now come back to the general B case and rewrite explictly the cavity method prediction for the free-energy per 
site, that is obtained from Eq. (|25[) by substituting /Uj_»j = r\: 



f3f = ln(5>(«x)P>>(<x') 






In ]T r^^V™' 



(29) 



and for the probability law on one site and two adjacent sites 



±-w(cr) (5>(«r')e' w ] 



(+i 



■2- edge 



(30) 



Note that the expression (|29[) is variational, in the sense that its derivative with respect to r\ vanishes on the solution 
of equation (f2"6")l . 

B. A convenient representation of rj(cr) 

We turn now to the problem of the determination of the probability law rj(cr) solution of Eq. f|26[) . As this is a law 
on the probability of N s Ising spins, its complete characterization should involve 2^ — 1 real numbers. Such a direct 
representation becomes very soon impossible to handle when N s grows, and in particular in the continuous time limit 
N s — > oo. We shall however see that an alternative representation allows to bypass this difficulty. First we define a 
probability law p(cr\h) on the configurations of a ring of spins cr by 



Pi<r\h) = -^ w(a) e^ h , Z(h) = ]T w(a) e^ h , 

2(h) ensuring the normalization of p(cr\h). These definitions allow to rewrite Eq. ([26"]) as 



(31) 



r)(cr) = ^2 Vifi) ■■■v( (T i) p(<r\<ri H ^ °i) 



Z(*! + ■ ■ ■ + <Tl) 



~-l 



(32) 



Suppose now that one has an estimation of rj(er) given by a representative weighted sample of Atraj elements {cj}, 
that is 



Mraj 

T](cr) = 22 a * ^(°" _£r *) 
»=i 



(33) 
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FIG. 5: Phase diagram: the critical temperature T c as a function of the transverse field B. 



such that the weights a% add up to one. A new estimation of r\ (i.e. a new set of configurations ct\ and weights a[) can 
then be obtained by plugging this estimation in the left-hand-side of Eq. (f3"2)l . which leads to the following procedure. 
To generate each of the new Atraj representants of i], one repeats in an independent way the steps: 

• draw independently I integers ii, . . . ,ii in [1, Atraj] with probability a^ 

• set h = cTi x + ■ ■ ■ + cr^ 

• generate a configuration &[ according to the law p(a\h), and set a[ = Z(h) 

Once the Atraj new elements have been generated one just has to multiply the weights a[ by a global normalization 
factor, and the new estimates can be again plugged in the right-hand-side of Eq. (|5^|) to approach its fixed point 
solution. 

At this point the continuous time limit N s — > oo can be taken without any difficulty of principle: as long as (3 is 
finite the spin trajectories have only a finite number of discontinuities, and can then be easily encoded with a finite 
set of reals corresponding to the time they change values. Moreover one can easily show that in this limit drawing a 
spin trajectory from the law p(cr\h) corresponds exactly to the procedure we defined in Sec. Ill Bl fwe show explicitely 
this correspondence for the normalization Z(h) in App. IB"]) . 

This representation of the probability law i](cr) by a sample of representative elements is a widespread technique in 
the field of disordered systems^. We warn however the reader accustomed to the classical cavity method that we did 
not use it exactly in the usual way, as the classical equivalent of rj{cr) is a single real number, not a probability law. 

Before closing this section let us discuss the computation of the physical observables in the continuous limit, taking 
as representative examples the longitudinal and transverse magnetizations. These can be obtained by taking the 
N s — > oo limit of the formulas © , (|10j) , which yields 



1 N 1 f 3 

= lim - Y(af) = Y m(<t)4 / di a 



(t) 



N 



N 



li m — V(af) =y^ /i ( C r)— jfer) , 



(34) 



where in the right-hand-side /i(c) is to be computed from Eq. (|30p . the sums over a are understood as sums over 
continuous time trajectories, and we have defined j{cr) as the number of discontinuities of a(t) on the interval [0, 0]. 



Continuous-time results 



We present now numerical results for 1 = 2; the behavior is qualitatively the same for any I. We followed the method 
of resolution presented above, using A/traj = 10 trajectories. For a fixed temperature, we initialized the population at 
S = 0; then each trajectory is constant and its value is chosen in such a way that the average magnetization is equal 
to the classical magnetization that can be obtained from the solution of Eq. (f28|) . In this way we select one of the two 
possible ferromagnetic states, that we can follow by increasing gradually the transverse field. We increased B with a 
step dB = 10~ 2 and for each step we let the population equilibrate for 10 2 iterations and averaged the observables 
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FIG. 6: Free energy as a function of B for T = 0.1. The curves for T — 0.3,0.2,0.1 coincide within our numerical precision. 
Therefore we assume that this curve is representative of the ground state energy as a function of B. The asymptotic values for 
B — > (ecs = 3/2) and for B — > oo (ecs = — B) are reported as dashed lines. The subleading correction for large B is oc _B _1 
and is not reported. 



over 10 3 subsequent steps. We checked that the weights a t defined in (|33p remain quite uniformly distributed over 
the population at all investigated temperatures and transverse fields. 

In figure [4] we plot the magnetizations m z and m x as functions of the transverse field B for three different temper- 
atures T < T C (B — 0) = l/arctanh(l/2) = 1.820. . .. At these low temperatures, the system undergoes a continuous 
phase transition at a critical value B C (T) of the transverse field. The transition is characterized by the vanishing of 
m z and by a jump in the derivative of m x . Unfortunately obtaining reliable values of m z close to B c is a difficult 
numerical task due to critical fluctuations that cause strong finite size effects in the size of the population TVtraj- 
Therefore we located the transition by the jump in the derivative of m x ; we fitted m x by linear laws close to B c from 
above and below, and determined their intersection. 

The resulting phase diagram is reported in fig. [5] We found that for T < 0.3 the temperature dependence of all 
observables is very weak and B c ~ 2.232, that is a reliable estimate for the T — ► limit and is in excellent agreement 
with the value determined in2^. The scaling of B c for T — > is compatible with the essential singularity that is found 
in the Curie- Weiss model, see Appendix [A] and in particular Eq. (jA8|) . The very weak (practically unobservable) 
temperature dependence of the free energy below T = 0.3 make us confident that the entropic term is negligible at 
these low temperatures and the free energy for T = 0.1, plotted in figure^is representative of the ground state energy 
eas- The latter has a weak singularity (discontinuity of the second derivative) at B c which is not observable with our 
numerical precision, but is easily observed by a direct computation of m x = — fs^i see figure 0] 

Overall, our results are in very good quantitative agreement with the ones obtained at T = ir>2i by a matrix 
product state description, except for the value of the exponent (3 characterizing the vanishing of the magnetization 
close to B c , m z ~ (B c — B)@ . Even if we do not have very precise data close to B c due to finite population-size effects, 
our data are compatible with the mean-field exponent f3 = 0.5 at all investigated temperatures, while in^l a slightly 
smaller value of /3 ~ 0.41 is reported at T = 0. We will further comment on this discrepancy in section HV CI 



D. Comparison with approximate treatments 



In this section we compare the previous results with approximated solutions of Eq. (126 
N s case, then we consider variational approximations to the solution. 



First we consider the finite 



Resolution at finite N s 



At finite N s , the distribution rj(cr) can be encoded with 2^" — 1 real numbers (e.g. the probabilities of each of the 
2 Ns configuration, with their sum constrained to be 1). Then Eq. (I26p can be rewritten as a fixed point equation for 
these numbers, and the solution can be found by simple iteration (below we refer to this procedure as "exact solution" 
for finite N s ). This method is extremely precise but computationally very heavy, as the time to solve the equation 
scales exponentially in N s . We are then limited to N s < 13; the computation for the largest N s = 13 took two weeks 
on a standard workstation, while we estimated the one for N s — 14 to take at least two months. 

Therefore, in order to study the approach to the limit N s — > oo, we resorted to the population method described in 
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Nb ) and solved Eq. 
oo. Similarly to the 



Sec. IIIIB1 We represented rj(cr) by a population of discrete time trajectories a = (a 1 , ■ • • , a 
by iteration following the procedure detailed in section IIIIBI without taking the limit 7V S 
continuous time case, we used A/traj = 10 5 and performed 10 2 iterations to achieve convergence, after which data have 
been collected along 10 3 iterations. The computation time is now polynomial in N s so we can go to much larger values 
(at the price of numerical precision due to finite A/traj)- Note that the continuous time computation is much more 
efficient. In the discrete time case, at low B and large T, the spins a a are constant along large time intervals and 
the information encoded in a discrete trajectory is redundant. On the contrary, at low T and large B the method at 
finite N s is obviously incorrect because of the discretization. One could think to adapt N s in order to find the optimal 
value for a given T,B; however this is most naturally done in the continuous-time framework where the number of 
spin flips is the natural variable describing a trajectory. Therefore we think that the finite N s population method is 
useful only for the illustrative purposes of this section. 

Data reported in this section and in figures [5] have been obtained by exact solution for N s < 13 and by the 
population method for N s > 13. We do not report detailed results for the magnetization as a function of B; as 
expected, the agreement with the continuous time computation is very good at high temperature, and becomes poorer 
and poorer on lowering the temperature. We focused on two points in the phase diagram; the first at intermediate 
temperature T = 1.0 and B — 1.5; the second at very low temperature T = 0.1 and B — 1.8. In figure [JJ we plot 
the longitudinal magnetization m z for different values of N s . As expected, finite N s effects are much stronger at low 
temperature. Note that for large iV s the leading correction is oc N s ~ 2 , as in the Curie-Weiss model (see Appendix [A"|l . 
For T = 1 deviations from the leading term are very small and N s ~ 10 is enough to get a fair estimate of the true 
m z . On the contrary, for T = 0.1 deviations are very large and N s > 256 is needed. 

In figure[S]we report the critical temperature as a function of B for different values of N s together with the continuous 
time result. Note that in the limit of large transverse field, each of the N s time slices becomes independent of the others. 
The model reduces to iV s copies of the classical system (B = 0) with a ferromagnetic coupling rescaled by 1/N S . Then 
in this limit the critical temperature is given by the classical one divided by N s , T C (N S , B — > oo) = l/[7V s arctanh(l/7)]. 
Below this temperature the system is always in the ferromagnetic phase, so that the quantum phase transition cannot 
be studied within this approximation. 

For generic spin-1/2 models, the finite N s approximation might be useful in order to understand the behavior 
at intermediate temperatures, and might also give quantitatively accurate results far from T = 0. This might be 
useful in cases where more complicated cavity treatments (such as the so-called one-step replica symmetry breaking) 
are necessary. Still the continuous time method appears to us preferable as it remains reliable down to very low 
temperatures. 



2. The static approximation 



A different strategy to obtain approximate solutions of these models is to use the variational property of the Bethe 
free energy (f2"9"]l : the free energy of the model is the minimum of the free energy defined by Eq. (f2"9")l over the function 
t]((t) (in the classical case the correctness of this procedure has been proven ini 3 -). Then one can propose a variational 
form for rj(cr) and minimize the free energy with respect to the free parameters to obtain an upper bound to the true 
free energy of the problem. 
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FIG. 8: Phase diagram at finite N a : the critical temperature T c as a function of the transverse field B. Note that for a fixed 
TVs, T C (B — » oo) = T C (B = 0)/N a . Filled squares are the result of the continuous time method, already reported in figure [5] 

A popular variational form is the so-called static approximation, that in our context amounts to the following ansatz 
for 77(c): 



77(0-) = n 



1 f 13 



'/ 






(35) 



i.e. the probability of a trajectory depends only on its average spin value along the imaginary time. This approximation 
has been widely used in computations based on the replica method for fully connected model a 14 i 15 ' 16 i 18 and has been 
recently applied to the random fc-satisfiability in a transverse fields. The present derivation based on the cavity 
method gives back the same results originally derived ir>22, but is simpler because the use of replicas is avoided. Here 
we discuss only the case of the ferromagnet on a regular graph but the equations can be easily generalized to more 
complicated cases where fluctuations of the couplings and/or the connectivity are present^. 
Equivalently we can rewrite the equation above as 



77(0-) 



N s 



dhp(h) ]J 



2 cosh h 



N e 



dmp(m) I 



mar. 



(36) 



where with a slight abuse of notation we used p for the distribution of both the effective field h and its associated 
magnetization m = tanh h. This expression makes evident that the assumption of the static approximation is that 
the spins in a trajectory are uncorrelated along the imaginary time, and subject to a common field extracted from 
the distribution p(h). 

From Eq. (|3"6"|) it follows that, in the limit N s — ► 00: 
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where in the last line we defined the function 



«'s 



E« 



(38) 



14 



m. 



1. 

0.8 






ttuL. 


i i 


i i 


4 

3.5 

3 


0.6" 


■■* 


^_ 




T \ 


^E_ 


2.5 


0.4 
0.2 




\ 


t 


\* \ 
\* \ 

\* \ 

i U ! 


M7 
AV 

4 V 
A V 

v - 

A 
V 


p(m) 2 
1.5 

1 

0.5 





0.5 


1 


1.5 2 
5 


2.5 " IP 




FIG. 9: (Left panel) Longitudinal magnetization m z as a function of the transverse field B. The full lines are the exact results 
reported in figure [4] for T — 1.7,1.0,0.1 (from left to right). Points are the results of the static approximation: T = 1.7 
(squares), T = 1.0 (circles), T = 0.1 (triangles), T = (open triangles). (Right panel) The function p(m) defined in the 
text at temperature T — 0.1 and transverse field B — 3.2,2.2, 1.2 (from left to right). For smaller values of B, the maximum 
approaches m = 1 and the weight of the S(m — 1) starts to increase until p(m) vanishes at B = 0. 



The following explicit expression for w s (m) can be found in (where it was called e 



-f3u {m)\ 



w s (m) = 



(3B 



VT 



- T h((3B^l-m 2 ) + S(m - 1) + S(m + 1) 



(39) 



with I\ the modified Bessel function of the first kind. For completeness we provide a proof of this result in App. [C] 
From the equations above one obtains the free energy (|29[) as a functional of p(m), that has then to be minimized. 
Remarkably, the resulting expression for the free energy is equivalent to the cavity free energy for a classical system 
whose variables are the continuous TOj G [—1,1] and whose Gibbs measure is defined by 






A ! 



Y[w s (mi) 



(40) 



i.e. it is obtained from the quantum measure ([5]) by replacing the quantum operators by their average and the 
transverse field term by its average as defined in (|38|) . The analogy with a classical system, or a direct differentiation 
of the free energy with respect to p(m), allows to write a cavity equation for p(m) which is the analog of (|26p : 



p( m ) = - Ws (m) ( f dm' p(m')e Pmm ' J 



(41) 



Given the structure of w s {m) it turns out that p{m) can be decomposed as a+6(m — 1) + a,-5(m + 1) + p(m), where 
p has its support strictly between —1 and 1. We solved Eq. (|4"Tj) by iteration, using a discretized representation 
of the regular part p and keeping explicitly the weights a±. Note that at B = the first term in (J39)) vanishes. 
Then a± = (1 ± tanh(/3/i))/2 where h is the solution of Eq. (|28|) : the static approximation is obviously exact in the 
classical case. For a given temperature T < T C (B — 0), we initialised p(m) on the classical solution and increased 



We can thus use the 



gradually the transverse field B, at each step iterating Eq. (|4ip until convergence (typically after ~ 10 iterations 
with a discretization step dm ~ 10~ 3 ). 

Eq. (|4i"|) can be further simplified in the limit T — > 0. As I\{x] ~ e 21 for large x, the Dirac distributions can be 
neglected in (|39f for any _B > 0, and at the leading exponential order w s (m) 



simplified expressions p(m) ~ e 



/3C(r, 



^i 



s /3BvT: 
; _,3e! at this order, and hence obtain from (HI 



C(m) = ei + By 1 — m 2 + I max[mTO' + C( m ')] 



(42) 



where the normalization ei must be determined in such a way that max m £(m) = 0. Like Eq. (|41|) . the latter equation 
can be solved iteratively. In this case a good starting guess is (,o(m) = Byl — m 2 + h$m + const, where the symmetry 
m — ► —m is initially broken by the term horn (with hg reasonably small) in order to select one state. Also in this case 
we used a discretization step dm ~ 10 -3 and observed convergence after ~ 10 2 iterations. 

The results of the static approximation are reported in figure |H1 In the left panel we compared m z obtained from 
the static approximation with the exact one obtained by the continuous time solution of the cavity equation. As 
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FIG. 10: Quantum heat bath (full symbols) vs. Loop Algorithm (open symbols) for a random regular graph of degree 3 with 
10 J spins at B = 1.6 (circles), in the ferromagnetic phase, and B — 2.6 (squares), in the quantum paramagnetic phase, for 
T — 1. Dashed lines correspond to the values computed with the cavity method. 

expected the approximation is very good for T > 1. and becomes poor close to T = 0. Still, the qualitative prediction 
of the static approximation remain reliable down to T = 0, even if the value of Bl tatlc (T = 0) = 3.0 predicted by 
the static approximation is different from the exact one, B C (T = 0) = 2.232. In the right panel of figure [9] we show 
the typical shape of p(m): it is a symmetric function at large B, where a± = 0. On lowering B below Static ; the 
distribution becomes asymmetric and a + > a_ > 0, even if a± remain very small at intermediate B. On approaching 
B = 0, a± grow faster and p(m) vanishes until the classical solution is recovered. 

In summary, the static approximation is a reliable variational tool to study qualitatively the phase diagram of the 
system. Remarkably, it can be solved down to T = (the solution at T — being easier than for finite T). 

IV. QUANTUM MONTE CARLO SIMULATIONS 
A. Methods and algorithms 

The quantum Monte Carlo method is an important tool that is widely used to study quantum statistical physics 
and quantum phase transitions, especially in the context of lattice models^. In this section, we discuss the application 
of the ideas discussed in this paper to quantum Monte Carlo simulations. 

1. A generic quantum heat bath Monte Carlo scheme 



The procedure we discussed in sec. Ill Bl to generate the "continuous time" spin configurations can be actually used 
directly as a continuous time quantum heat-bath algorithm. Once the procedure that generates the new continuous 
time configuration given the local field trajectory is available, the implementation of the Monte Carlo simulation is 
rather straightforward: one just randomly picks a site, computes the local field trajectory due to its neighbouring 
spins, and generates a new imaginary-time trajectory of this spin according to the rules discussed in sec. Ill Bl 

The advantage of this method is that one can apply it to any discrete spin models on any kind of lattice. In the case 
of ferromagnetic non frustrated interactions, there exist of course a number of algorithms that allow to considerably 
speed up the simulation and this is the subject of the next section. However, for highly disordered and frustrated 
systems such as spin glasses, for instance, no such generic cluster algorithm exists even at the classical level (although 
interesting progresses are being made, se o 49 ' 50 i 51 ' 52 ) and most classical simulations^ rely on the Metropolis or the 
heat bath algorithm and the parallel tempering technique^. In the case of quantum spin glasses, to the best of our 
knowledge, only the usual discrete time Suzuki- Trotter decomposition has been trie d 55 i 56 ' 57 . Another important case 
where no loop algorithm is known is the one of multi-spin interactions, where the problems are defined on a factor 
graph (see section |V| , which is very common in the context of constraint satisfaction problems such as satisfiability 
or in coding theory There is thus a clear need of an efficient generic heat-bath strategy when no cluster or loop 
algorithm exists, that would result in a substantial gain of simulation time. 

In Fig. 1101 we show the first application of our continuous time quantum heat-bath algorithm to the case of the 
ferromagnetic model on a regular random graph of connectivity 3, as studied in the rest of the present paper. The 
results of our algorithm converge fastly, with respect to the number of Monte Carlo sweeps, to the asymptotic ones 
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FIG. 11: Longitudinal and transverse magnetization curves for temperature T = 1 (left) and T = 0.1 (right). Continuous time 
cavity method (m z , solid line and m x , dashed line) vs. Monte Carlo simulations for different sizes increasing from the top to 
the bottom (N — 64, 128, . . . , 2048). For the largest size the agreement with the cavity result is excellent (except close to the 
critical point). 

for very large samples (i.e. N = 10 5 spins). In the paramagnetic phase, the convergence is even faster than for the 
loop algorithm discussed in the next section. The application of these ideas and methods to more complex models is 
an interesting direction of study. 



2. Continuous time loop algorithm 



Since we are dealing with a ferromagnetic, non frustrated model, we expect the usual cluster and loop 
algorithm o 36 ' 39 ! 58 to allow for a significant speed up of the simulation with respect to the heat-bath procedure, 
even on random graphs. We have thus implemented the algorithms devised by Rieger and Kawashima in 3 ^, which is 
an adaptation of the classical Swendsen-Wang 59 algorithm to quantum systems in continuous time. The results are 
compared with the heat-bath method in Fig llOl Indeed, the loop algorithm performs very well and better than the 
heat-bath method close and below the critical point, and we thus have used this method to compare the results of 
the cavity approach with finite size instances. 



B. Cavity method versus simulation 

How does the cavity predictions compare with numerical simulations? To answer this question, we have performed 
quantum Monte Carlo simulations on random regular graphs of N spins and connectivity 3. For large N we expect 
the results to be self- averaging and thus we always consider a single instance, and do not perform averages of many 
realizations. The results for two different temperature T = 1 and T = 0.1 are shown in Fig llll where we plot the 
longitudinal and transverse magnetization as function of B for different sizes from N = 64 to N — 2048. The 
agreement with the asymptotic cavity result is perfect (apart from finite size effects in m z for B s=s B c , see next 
paragraph): this demonstrates the correctness of the approach we have developed in the present paper. 



C. Finite size scaling and critical exponent 



The transition between the ferromagnetic and the paramagnetic phase is of second order; below the threshold B C (T) 
the longitudinal magnetization grows with an exponent /3, i.e. m z oc (B c — B) , in the thermodynamic limit. We 
have used finite size scaling technique o 60 i 61 to analyze our data in the neighborhood of the transition and to check 
the mean-field value of the exponent (3 — 1/2. Let us first briefly recall the basic idea of the finite size scaling method 
and the way in which it has to be amended for mean-field models. 

In a generic infinite size d-dimensional model, in the vicinity of a second order phase transition driven by a parameter 
denoted B, the correlation length diverges as £ ex \B — B c \~". For a system of finite extent L the observables depends 
on the size through a scaling function of the ratio L/£. This has to be corrected for dimensions d larger than the upper 
critical dimension d n of the considered universality class, or when the model lacks any underlying finite dinrensional 
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structure. In that case the scaling function is foun d 62 ' 63 ' 64 ! 65 ' 66 to depend on the size N (total number of degrees of 
freedom, equivalent to L d of a d-dimensional model) through JV 1 /^"^ (B — B c ), where v takes its mean-field value in 
the universality class under investigation. More explicitly the scaling forms of the longitudinal magnetization and of 

the Binder cumulant g = 



<mi>* 



read 



i z (B,N) = N-P'^m (V/(^) (B - B c j\ , 
g(B,N) = g (jyVC*.") (B - B c j) ■ 



(43) 



We present in Fig. [T2] the results of such an analysis for our Monte-Carlo data obtained at temperature T = 1. At this 
positive temperature the critical behavior of a quantum d-dimensional system is equivalent to a classical d-dimensional 
Ising mode l 13 ' 29 , which implies d u = 4 and v = 1/2. Using the cavity method predictions B C (T = 1) = 1.98 and 
(3 = 1/2 we found a very good data collapse (see Fig. [l"2"|) , which confirms the validity of these values of (3 and B c . 

To close this section, let us discuss the value of the critical exponent j3 at T = 0. We believe that because of 
the mean-field nature of the Bethe lattice model, (3 keeps the same value 1/2 at positive and zero temperatures. 
This is indeed what happens for the Curie- Weiss model, see Appendix [AJ A simple argument is the following. For 
ferromagnetic models the Suzuki- Trotter formalism and numerical simulations suggest that the critical behaviour 
of the quantum d-dimensional model at zero temperature corresponds to the one of the classical model in d + 1 
dimension s 13 ' 29 ^ 8 . Hence if a model behaves in a mean-field way at positive temperature it should also do so at zero 
temperature, having formally gained one more spatial dimension. To give further credit to our theses we performed 
quantum Monte-Carlo simulations at very small temperatures, T — 0.01. The plot of Fig. IT51 shows again a very good 
collapse with the value j3 — 1/2, using the zero temperature critical value of the transverse field extracted from the 
cavity computation, B c = 2.232. This collapse has been obtained with d u = 3; indeed at such a low temperature we 
are well inside the low temperature regime where the transition is truly quantum, hence the shift of the upper critical 
dimension to account for the imaginary time supplementary dimension. 

A different value of the zero temperature exponent, f3 = 0.41, has been obtained inSi using the matrix product 
state ansatz for the description of the groundstate. Though we cannot rule out the possibility that a crossover occurs 
at temperatures even lower than T = 0.01, we believe that (3 — 1/2 down to T = 0, in accordance with the mean-field 
nature of the Bethe lattice, and that the different value reported in 2 ^ might be due to the truncation of the matrix 
product state ansatz. 



V. THE GENERIC REPLICA SYMMETRIC QUANTUM CAVITY METHOD 



In this final Section we give a more generic description of the replica symmetric quantum cavity method. As the 
main ideas should have already been conveyed by the example of the ferromagnet on the random regular graphs we 
shall mainly emphasize the differences and complications that arise in more general models. 

Let us consider the class of Ising spin models whose classical energy E(<r) can be decomposed into the sum of M 
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FIG. 13: Finite size scaling analysis of the longitudinal magnetization at T — 0.01, for sizes N = 32, 64, . . . , 2048. 



interaction terms a = 1, . . . , M, each of them depending on a finite number of spins, 

M 



E{a)= y £ d s(a Sa ,J a ) 



(44) 



In this expression da is the subset of spin indices the a-th interaction effectively depends on, g_g a = ((Tj : i € da) is a 
shorthand for the configuration of those spins, and J a denotes coupling constants that might appear in the definition 
of the a-th interaction term. This decomposition is conveniently represented as a factor grapbM, i.e. a graph with 
two kind of vertices (see Fig. [14] for an illustration): squares stand for the interactions a — 1, . . . , M, while circles 
represent the variables i = 1, . . . , JV. An edge is drawn between an interaction a and a variable i whenever a depends 
on i, in other words whenever i G da. In this graphical representation da is thus the set of nodes adjacent to a. 
Similarly we shall denote di the set of neighbors of i, that is all the interactions that depend on oi. The extended 
Ising model defined in Eq. (O preserves the topology of the interactions of the classical energy E(a), the latter being 
reproduced identically in the various imaginary time slices, 

M 



E{a) = ^e{a dal J a ) , 

a=l 

1 ^ 



(45) 



while the transverse weights w[<Ti) are local in the spin indices i. 

As we did on the ferromagnet example we first state the solution of this extended model on a tree. The difference is 
that we have now two kind of "messages", from variables to interactions and vice versa. Let us denote /i a _>i(£Ti) the 
probability law of <Xj in the model corresponding to the factor graph where all interactions but a have been removed 
in the neighborhood of i, and similarly /Xj_» a (<T.j) for the factor graph with only a removed from di. These quantities 
obey the following recursion equations, 



IM-> a (<ri) = w(a-i) Yi fJ-b^i(cri) 



(46) 



b£di\a 



ia->i(tn) = Y] 



,-Pz(SL3a> J o.) 



n ^^K") 



{"j}j£ 



j£da\i 



on all edges of the factor graph, the various constants z being normalization factors. The reader will easily verify that 
these equations reduce to Eq. I|22|l in the ferromagnetic case with e(ai,<Tj) — —cricrf, in this case we could eliminate 
one type of message (from interactions to variables), as the interactions were only pairwise. 

We shall now consider ensembles of random factor graphs, denoting E[-] the expectation over the distribution of 
the graphs and coupling constants. We assume all interaction nodes to involve a fixed number k of variables (the 
ferromagnet corresponded to k = 2), while the degree distribution of the variables is specified by a probability law 
qd over the positive integers (all random graphs verifying these constraints are equiprobable in the ensemble; the lack 



19 



> 



# 



FIG. 14: An example of a factor graph 



of a finite dimensional a priori structure is the origin of the mean- field character of this family of models). We shall 
denote "fk = ^2 d dqd the average variable degree. Moreover the M = jN coupling constants J a are drawn in an 
identical, independent way for each of the interactions. Suppose the recursion equations (|46|) are solved on a factor 
graph sampled at random from the ensemble under consideration, and that an edge from a variable to an interaction, 
call it i — ► a, is chosen uniformly at random. The probability law /Uj_> it bears is itself a random variable, denoted 
rj in the following, due to the random choices of the graph and of its coupling constants. This random variable can 
be related to its neighboring equivalents by the local relations (|46p . Suppose that i has d adjacent interactions apart 
from a (see Fig. [15]for an illustration). Then one has, from (|46|) . 



r/(cr) — —w(cr) 



{<T a 




i{Ca,i) 



-PY. d a = 1 e{<T,<r a 



t-uJa) 



(47) 



where the r\ a ^ are the d(k — 1) probability laws that determine r\. The hypothesis of the replica symmetric cavity 
method is to assume that the above form is correct in spite of the graph not being globally a tree, and that the r\ a ^ 
are independent identically distributed copies of the random variable rj, which is written in formula as 



n = f(v 



i,i) 



>»7i,fc-i> 



)%,1, ■ ■ ■ ,Vd,k-l)Jl) ■ ■ -,Jd) , 



(48) 



the symbol = denoting identity in distribution of random variables, and the function / being an abbreviation of the 
right-hand-side of (f47|) . Note that the connectivity random variable d is not distributed according to qd but rather 
Qd = [d + V)qd+i / lk: as one uniformly choose a random edge of the graph, and not a random site, this sampling 
favours larger connectivity variables. For the ferromagnet on the random regular graph one had qd = 8d,i+i while 
qd = 5d,i- 

The physical observables of the system can be computed from the solution of this distributional equation. The 
thermodynamic limit of the free-energy per site is found to be 



(if = lim — EflnZl = E 



r} a ,i{ (T a,,i) w{cr)e 




-PT,i =1 £(° 



-iJa 



7(fc-l)E 



-/3e(<Ti,...,<T fc ,J) 



(49) 



where the expectations of the left-hand-side are over the choice of the random factor graphs, and for the right-hand- 
side over independent copies of the random variable r\ and the coupling constants J, while in the first term d is drawn 
from the law qd- This is a generalization of the expression l|29|) found for the ferromagnet. Similarly the average 
marginal of the law ([5]) for an arbitrary site reads 



fi(cr) = E 



— w(cr) 



E 

, ,ie[i,k-i] 



Q%,i(o"o,») 



,-/3Eo=l e{<r,<ra,u-><ra,h-i,Ja) 



(50) 
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FIG. 15: A pictorial representation of Eq. (|47|l . 



with d drawn from qd, while the average marginal law of the k variables in an arbitrary interaction is 



/i(cri,...,<Tfe) =E 



i(n *(»*))«- 



■/9e(«Ti 



,Tk,J) 



(51) 



the factors 1/z in these last two equations ensuring their normalization (cf. Eq. (|30[) for the regular ferromagnet). 

Let us now discuss a method of resolution of the distributional equation (|48|) . consisting in encoding the distribution 
of 77 by a sample (or population) of a large number J\f of representatives 77^ . From an arbitrarily initialized population, 
one applies iteration steps according to Eq. (|48|) : one draws an integer d from q^, d coupling constants J\,...,Jd 
and d{k — 1) indices in [1,A/]. The right-hand-side of Eq. (|15|) is then computed from the corresponding d(k — 1) 
representants of r\ randomly chosen in the population, and the resulting 77 is used to replace one discarded element of 
the population. The iteration of these steps brings the population close to a fixed point of the distributional equation; 
then physical observables like the average free-energy (see Eq. (|49)l ) can be computed, evaluating the expectations 
over the random variable r\ as a sampling from the approximate representation provided by the finite population. 
Compared with the classical replica symmetric cavity method the difficulty is that each of the r\{ in this population is 
itself a probability distribution over configurations of spin rings, or in the continuous limit over spin trajectories. It 
is however possible to apply the trick explained in the simpler case of the ferromagnet, and encode each of the rjt as 
a sample of A/traj spin trajectories. Let us introduce some further notations in order to give an explicit form of the 
updating procedure of the population of r/j's. 

The dependence of e(-, J) on one of the Ising spins can always be parameterized through two functions u and v as 



e(a,ai,...,a k -i, J) = - au(ax,...,a k -i,J)+v(ai,...,ak-i,J), 
u playing the role of an effective magnetic field acting on a. For JV s -fold replicated variables we define 



(52) 



tt(<Ti, 



,<T k -i,J) = (u(al,...,a k _ x , J),...,u(CTf a ,...,crf_! 1 , J) , 

,0-fe-l,,/) = — $3uK , ,...,<7fc_i, J) 



a=l 



Using the definitions in Eq. (f3i~j) . we can then rewrite Eq. (|47|) as 



(53) 

(54) 



1]((T) 



e (n »?«.<(*«.<)) P {*\h{{v a „j a })) z{{(Ti >f a}) 



with 



h({(Ja„i,Ja\) = ^2u((T a}1 ,...,er atk _i,J a ) , z({er a4 ,J a }) = Z(h({(T a , t ,J a }))e 



-0^2a = l v i cr a,,l, — ,Va,k-l,Ja 



(55) 
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We can thus apply the sampling procedure explained in Sec. IIIIB1 the only change being the different definition of 
the effective field trajectory h, and the inclusion of a contribution arising from the function v in the weights of the 
generated spin trajectories. Let us emphasize that at variance with the classical cavity method, here we have to deal 
with a population of population of spin trajectories (of total size Af x A/t ra j) already at the replica symmetric level. 

VI. CONCLUSIONS 

The explicit procedure for the construction of continuous time spin trajectories presented in Sec. Ill B I allowed us to 
make progresses both on the analytical side, with an improvement over the discretized time cavity method of 21 and the 
static approximation o^, and on the numerical simulations side, with a generic quantum Monte Carlo procedure for 
spin-1/2 models in transverse field. The ferromagnetic model we studied in this article is only the simplest of a large 
family that can be tackled, at the price of an heavier computational cost, with the same methods. We believe that its 
study was in any case worthwhile from a methodological point of view: its simplicity permitted a complete resolution 
of the quantum cavity equations, which (i) showed a perfect agreement with Monte Carlo simulations, hence giving 
credit to the conjecture that the replica symmetric cavity method leads to exact results in the thermodynamic limit 
for sparse mean-field ferromagnetic systems, as was proved for classical models ini£; (ii) allowed to test quantitatively 
some approximate treatments (finite N s and static approximation). 

The more interesting cases we plan to address in the future by mean of the cavity method will involve fluctuating 
connectivities, in order to study the role of these local fluctuations on the critical behaviour of the models, glassy 
phases at low temperatures (a particularly motivating case will be the regular multispin ferromagnet studied at the 
classical level in 4 ^), and models related to quantum computing issues, as the random fc-satisfiability model in a 
transverse field22. Note that already the replica symmetric treatment of these models will involve a population of 
populations of spin trajectories, as explained in Sec. |Vl which will make an exact treatment of the one step replica 
symmetry breaking version of the cavity method extremely challenging. We hope however that the better control 
of the approximative treatments we gained on the simple ferromagnet will help to devise appropriate approximate 
strategies to handle these cases. 

At the same time, we proposed a heat bath Monte Carlo strategy that should be applicable to general spin-1/2 
models. The performance of this algorithm should be tested on frustrated models for which cluster algorithms are 
not easy to implement, and we expect that in such cases the heat bath method could give interesting results. 

The strategy developed above become inefficient at very low temperatures, because the spins jump many times 
along a trajectory and the amount of information needed to encode t](<t) by a population of trajectories becomes very 
large. Therefore an important directions of research would be to look for possible simplifications of the formalism 
in the limit (3 — » oo. Moreover, one could use the recursion equations l|3o]) on a given sample, thus constructing a 
quantum belief propagation algorithm 6 -^. Finally, it is worth to note that the procedure described in section III Bl on 
which our results are based, is in principle not restricted to spin-1/2 models. It should be possible to generalize it to 
any model built on discrete degrees of freedom^, for instance hard-core bosons or higher spin models. 
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APPENDIX A: QUANTUM CURIE- WEISS MODEL 

This Appendix is devoted to a study of the simplest quantum mean-field ferromagnet, namely the fully-connected 
Curie- Weiss quantum mode l 44 ' 45 ! 46 . The model is defined by the Hamiltonian 

j N N N 

£=-2iv£<* , i- B 5>?- fc £ ?' (A1) 

i,j=l i—1 i— 1 

where the scaling of the coupling constant is chosen appropriately to make the thermodynamic limit well-defined. 
At variance with the Bethe lattice model here each spin interacts with all others. Applying the Suzuki- Trotter 
decomposition described in Sec. Ill Al leads for a finite number JV" S of Suzuki- Trotter slices to the following expression 
of the partition function: 
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FIG. 16: Left: the phase diagram of the quantum Curie- Weiss model. Right: longitudinal and transverse magnetizations as a 
function of the transverse field, for T = 0.7. 



A' 



z = e ( n u, (° - *) exp 



i=l 



0h 

iVs 



AT S 



E-f 



exp 



(3J 

2NN, 



N s N 



a — 1 i,i=l 



(A2) 



Wc then perform N s Hubbard-Stratanovitch transformations to disentangle the quadratic terms and obtain 



/ j3JN \ 



* „ n s 



yn- 



m exp 



-TV 






N e 



Y,(m a ) 2 + NlnTr [l[e^ h+Jma ^ e^ B °* 



(A3) 



Evaluating these integrals by the saddle-point method in the thermodynamic limit and selecting the cyclically invariant 
saddle point yields 



lim — In Z = sup 

AT^oo N m 



0J 2 

rr 

2 



(0 



m + In ir e "= v ' e N = 



A. 



(A4) 



This can be further simplified if the N s — ► oo limit is performed afterwards and yields for the free-energy per site 



/ = inf 



J 



— — In ( 2 cosh 
P 



^(h + Jmf + B 2 ^ 



(A5) 



Using the variational character of this expression the longitudinal and transverse magnetizations (i.e. m z and m x ) 
can be obtained by taking the explicit derivatives with respect to h and B respectively, 



h + Jm 



m x = o\ 



y/(h + Jm) 2 + B 2 

B 
y/(h + Jm) 2 +B 5 



tanh (Pyj{h + Jm) 2 + B 2 \ 
(l3^/(h + Jm) 2 + B 2 ^ 



tanh 



(A6) 

(A7) 



where m is taken as the solution of the saddle-point equation. The latter is easily found to imply that m = m z at the 
saddle-point. In the following we take J = 1 to simplify the notations. 

In absence of the transverse field (B = 0) one recovers the classical Curie- Weiss model, with m — m c \{j3, h) solution 
of the traditional equation m = tanh(/3(m + h)). The ferromagnetic transition is signaled by the appearance of a 
non-trivial solution at h = 0, which is possible for small enough temperatures, i.e. for (3 > f3 c (B = 0) = 1. 

Consider now the solutions of the saddle-point equation with B > and h = 0. The paramagnetic solution m = 
exists for all temperatures and transverse fields. For a strictly positive solution m the saddle-point equation reduces 
to \Jm 2 + B 2 = tanh(/3\/m 2 + B 2 ), that is m(0,B,h = 0) = y/m c \{(i, h = 0) 2 — B 2 . This is possible only for small 
enough temperatures and transverse fields, such that the argument of the square root remains positive. The line of 
transition in the (B, T) plane is such that B c (j3) = m c i(/3, h = 0), see the top panel of Fig. [TBI Note the vertical slope 
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of the transition line in the neighborhood of the quantum critical point in (B = 1, T = 0); in fact one can perform an 
asymptotic expansion in this region, to show that 

B c {T) T ^l-2e- 2 P , (A8) 



or equivalently, 



T C {B) ~ — 7-i r . (A9) 



The bottom panel of Fig. [16] shows the evolution, as a function of the transverse field, of the longitudinal and 
transverse magnetizations at a temperature smaller than the classical critical temperature T C (B = 0) = 1. The 
transverse magnetization is continuous but its derivative has a finite jump at the transition. One can indeed show 
from Eq. (|A7[) that m x = B in the ferromagnetic phase, while m x = tanh(/3-B) for the paramagnetic one. 

Another thermodynamic quantity easily computed for the Curie- Weiss model is the ground-state energy per spin, 
e gs (_B), obtained from (|A5|) in the limit j3 — > oo. One finds 

^>={:f +B2 > tlif-l 

which shows that e gs (B) and its first derivative are continuous at the transition, while its second derivative has a 
finite jump. 

Let us finally argue about the scaling of the finite N s corrections, reconsidering the step we took between Eqs. (|A4[) 
and (|A5[) . At the next-to-leading order one would have obtained 



Tr(7e£ ( ' l+Jm)CTi e£ B<TX ) ArsN ) =TiYexp [3(h + Jm)a z + $Bo x + — f3 2 B(h + Jm)[a z , o x ] (l + 0(AT 2 ))") . 

One can then check explicitely that the eigenvalues of the matrix in the exponential are not modified at the order 
Ag -1 , hence the finite N s correction on the observables should be of order N s ~ 2 , as we observed in the study of the 
Bethe lattice ferromagnet. 

APPENDIX B: AN IDENTITY IN THE CONTINUOUS TIME LIMIT 

In this appendix we briefly justify the claim made in Sec. IIIIBI of the equivalence in the continuous time limit of 
the definition of Z(h) given in (|3Tj) with the one of (fl9|) . Before taking the N s — > oo limit, the former reads 

N s 

Z(h) = ^ ™(<r) exp[/3 <r ■ /i] = J^ n(<? a \e^ ha,7 *e& Ba *\(j a+1 ) . 

1 CT Wa a=l 



CT 1 , <7 



We shall denote N s ^> = N s \^ //3 the lengths of the constant longitudinal field intervals (see Fig. [2]), expressed in 
number of Suzuki- Trotter slices. Then 

z w = e nwt (,, )i(^ k " ia ^ B, f , "V(t (,+i) )) <7(i(p +i >) = <7(*<°>) 

<r(t(o)),...,o-(t(p))i=l 

E f[nt (l) )\e xi ^ h( ' ) ^ +B ^\a(t^+ 1 ^) . (Bl) 

<T(t(°)),...,cr(t<P))« = 



In the last step we have taken the A s — > oo limit to obtain the expression of Eq. (|T9 
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APPENDIX C: THE COMPUTATION OF w B (m) 



In this appendix we derive the expression (|39[) for the average transverse weight in the static approximation. 
Rewriting the sum over cr of (|38p as a sum over paths in the continuous time limit, one obtains 



I ^ vr*r* Fa* f u xl 2t 1 -2t 2 + ----2t n +0 \ 
,(m) = 2^/Z^ B dtl dt 2 ... dhn6[m-a - ), (CI) 



<5(m - 1) + 5{m + 1) + ^ ^(B/3) 2n / da* / dx 2 . . . dx 2n 



/3 

5 (2x\ — 2x 2 + • • • — 2x n + 1 — ma) 



where in the second line we have isolated the contribution of the constant trajectories, and changed variables from ti 
to Xi — U/ (3. Let us concentrate on the term corresponding to a given value of n > 0. We perform a further change 
of variables, setting yi — xi — Xi-\ (with xq = 0), the intervals between the reduced time of flips in the trajectory. 
The Jacobian of this change of variables being 1, the integral over x±, . . . , x 2n can be rewritten as 



/ dyi . . . dy 2n l(yi +y 2 -\ h y 2n < 1) S(2y 2 + 2y 4 -\ \- 2y 2n - (1 - ma)) 

Jo 

= i J dS dS c Pn (S )p n (S c ) 1(5 + S c < 1) S \S e - X —^° 

l + mo- 



(C2) 



77P™ ( s ) / dS p n (S) . 



Here we have used !(•) as the indicator function of an event, S = 2/1 + 2/3 H h Vin-\ and 5 C = j/2 + ■ • • + 2/2n the 

sum of the odd/even y^s, and p n (S) as the density of the distribution of the sum of n independent random variables 
uniformly distributed on [0, 1]. It is easy to prove by recurrence that for S G [0, 1] one has p n {S) = S 71 ^ 1 /(n — 1)!. 
Collecting these facts together leads to 

u, s = S(m -l) + S(m + l) + ± * 1 / B)2 "^:r 2) "' 1 • (C3) 

n=l '\ '' 

which is indeed equal to (|55|) thanks to the series expansion of the Bessel function. 
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